Classification of the pairing transition in finite fermionic systems 
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In this work we employ a simple pairing interaction model in order to study and classify an 
eventual pairing phase transition in finite fermionic systems. We show that systems with as few as 
~ 10—16 fermions can exhibit clear features reminiscent of a phase transition. To classify the nature 
of the transition we apply two different numerical methods, one based on standard thermodymanics, 
and another based on a recently proposed scheme by Borrmann et al. [Phys. Rev. Lett. 84, 3511 
(2000)]. The transition is shown to be of second order, in agreement with results for infinite fermionic 
systems. 
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The standard BCS theory has been widely used to de- 
scribe systems with pairing correlations and phase transi- 
tions to a superconducting phase for large systems, from 
the solid state to nuclear physics, with neutron stars as 
perhaps the largest object in the universe exhibiting su- 
perfluidity in its interior. An eventual superfluid phase 
in a neutron star will condition the neutrino emission and 
thereby the cooling history of such a star, in addition to 
inducing mechanisms such as sudden spin ups in the ro- 
tational period of the star; see, for example, Ref. jjj for a 
recent review. For an infinite system, such as a neutron 
star, the nature of the pairing phase transition is well 
established as second order. 

When a system of correlated fermions such as electrons 
or nucleons is sufficiently small, the fermionic spectrum 
becomes discrete. If the spacing approaches the size of 
the pairing gap, superconductivity is expected to break 
down however, recent experiments on superconduct- 
ing ultrasmall aluminum grains by Tinkham et al. || 
revealed the existence of a spectroscopic gap larger than 
the average electronic level density. This feature was in- 
terpreted as a reminiscence of superconductivity and re- 
newed the interest [Q-|7j in studies of what is the lower 
size limit for superconductivity. 

Other finite fermionic systems such as nuclei are ex- 
pected to exhibit a variety of interesting phase-transition 
like phenomena, like the disappearence of pairing at a 
critical temperature T c « 0.5 — 1 McV or the nuclear 
shape transitions of deformed nuclei associated with the 
melting of shell effects at T c « 1 — 12 McV. In recent the- 
oretical and experimental studies [^],[| of thermodynami- 
cal properties of finite nuclei, the heat capacity has been 
found to exhibit a non-vanishing bump at temperatures 
proportional to half the pairing gap. These bumps were 
interpreted as signs of the quenching of pair correlations, 
representing in turn features of the pairing transition for 
an infinitely large system. 



In the study of phase transitions in e.g., solid state, 
nuclear and high energy physics, it is important to know 
whether a given transition really is of first order, discon- 
tinuous, or if there is a continuous change in a physical 
quantity like the mean energy, as in phase transitions 
of second order. If one works in the canonical or grand 
canonical ensembles, for finite systems it is rather diffi- 
cult to decide on the order of the phase transition. This 
is due to the fact that in ensembles like the canonical, any 
anomaly is smeared over a temperature range of 1/N, N 
being the number of particles. In the analysis of finite 
systems, both a <5-f unction peak and a power law sin- 
gularity sharpen as the number of particles is increased, 
making it difficult to distinguish between the two cases, 
see, for example, Ref. j |To| . 

In addition, first order phase transitions in finite sys- 
tems have recently been inferred, theoretically and ex- 
perimentally, from observed negative heat capacities that 
are associated with anomalous convex intruders in the en- 
tropy versus energy curves, resulting in backbendings in 
the caloric curves; see, for example, Refs. Neg- 
ative heat capacities are often claimed to appear only 
in calculations done in the microcanonical ensemble and 
are thought to vanish in the canonical or grand-canonical 
ensembles. 

We aim in this Letter to identify the nature of the pair- 
ing transition, if any, in a finite fermionic system. Since 
we are dealing with pairing correlations, our Hamiltonian 
is 

H = ^ e t aja t — G T"' alalajaj , (1) 

i ij 

where and a are fermion creation and annihilation op- 
erators, respectively. The indices i and j run over the 
number of levels L, and the label T stands for a time- 
reversed state. The parameter G is the strength of the 
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pairing force while Si is the single-particle energy of level 
i. 

We assume that the single-particle levels are equidis- 
tant with a fixed spacing d. Moreover, in our simple 
model, the degeneracy of the single-particle levels is set 
to 2J+1 = 2, with J — 1/2 being the spin of the particle. 
Introducing the pair-creation operator — a\ m a\_ m , 
one can rewrite the Hamiltonian in Eq. (HI) as 

H = d^N l -GY / S+S-, (2) 

i ij 

where Ni = a\ai is the number operator. Seniority S 
is a good quantum number and the eigenvalue problem 
can be block-diagonalized in terms of different seniority 
values. Loosely speaking, the seniority quantum number 
S is equal to the number of unpaired particles. 

For systems with less than ~ 16 — 18 particles, this 
model can be diagonalized exactly, and we can obtain 
all eigenstates. In our studies below, we will always con- 
sider the case of half- filling, i.e., equally many particles 
and single-particle levels. This case has the largest di- 
mensionality: for 16 particles in 16 doubly degenerate 
single-particle shells, we have a total of 4 x 10 s states. 
We choose units MeV for the energy and set G = 0.2 
MeV in all calculations while we let d vary. 

Through diagonalization of the above Hamiltonian we 
can define exactly the density of states SIn(E) for an N- 
particle system with excitation energy E. An alternative 
to the exact diagonalization, would be to use Richard- 
son's well-known solution jll|, however, we are interested 
in all eigenstates, and the amount of numerical labor will 
most likely be similar. The density of states is an essen- 
tial ingredient in the evaluation of thermal averages and 
for the discussion of phase transitions in finite systems. 
For nuclei, experimental information on the density of 
states is expected to reveal important information on nu- 
clear shell structure, pair correlations and other correla- 
tion phenomena in the nuclconic motion. 

The density of states SIn(E) is the statistical weight 
of the given state with excitation energy E, and its log- 
arithm 

S N {E) = k B lnn N (E), (3) 

is the entropy (we set Boltzmann's constant fc^ = 1) 
of the TV-particle system. The density of states defines 
also the partition function in the microcanonical ensem- 
ble and can be used to compute the partition function Z 
of the canonical ensemble through 

Z(/3) = Y / ^N(E)e^ E , (4) 

E 

with P = 1/T the inverse temperature. With Z it 
is straightforward to generate other thermodynamical 
properties such as the mean energy (E) or the specific 
heat Cy. 



The density of states can also be used to define the free 
energy F{E) in the microcanonical ensemble at a fixed 
temperature T (actually an expectation value in this en- 
semble) , 

F{E) = -Tin [Q N (E)e-^ E ] . (5) 

Note that here we include only configurations at a par- 
ticular E. 

The above free energy was used by e.g., Lee and 
Kosterlitz Jig] , based on the histogram approach for 
studying phase transitions developed by Ferrenberg and 
Swendsen jlTj , in their studies of phase transitions of clas- 
sical spin systems. If a phase transition is present, a plot 
ofF(E) versus E will show two local minima which corre- 
spond to configurations that are characteristic of the high 
and low temperature phases. At the transition tempera- 
ture Tq the value of F(E) at the two minima equal, while 
at temperatures below Tc, the low-energy minimum is 
the absolute minimum. At temperatures above Tc, the 
high-energy minimum is the largest. If there is no phase 
transition, the system developes only one minimum for all 
temperatures. Since we are dealing with finite systems, 
we can study the development of the two minima as func- 
tion of the dimension of the system and thereby extract 
information about the nature of the phase transition. If 
we are dealing with a second order phase transition, the 
behavior of F(E) does not change dramatically as the 
size of the system increases. However, if the transition 
is first order, the difference in free energy, i.e., the dis- 
tance between the maximum and minimum values, will 
increase with increasing dimension. 

To elucidate the nature of the pairing transition we 
employ two different approaches, which both rely on our 
ability of computing the exact density of states fl(E). 

First, we calculate exactly the free energy F(E) of 
Eq. through diagonalization of the pairing Hamilto- 
nian of Eq. (|lj) for systems with up to 16 particles in 16 
doubly degenerate levels. For d/G = 0.5 and 16 single- 
particle levels, we develop two clear minima for the free 
energy. This is seen in Fig. [l] where we show the free 
energy as function of excitation energy using Eq. (^|) at 
temperatures T = 0.5, T = 0.85 and T = 1.0 MeV. The 
first minimum corresponds to the case where we break 
one pair. The second and third minima correspond to 
cases where two and three pairs are broken, respectively. 
When two pairs are broken, corresponding to seniority 
5 = 4, the free energy minimum is made up of contribu- 
tions from states with S = 0,2,4. These contributions 
serve to lower the free energy. Similarly, with three pairs 
broken we see a new free energy minimum which receives 
contributions from S = 0, 2, 4, 6. At higher excitation en- 
ergies, population inversion takes place, and our model 
is no longer realistic. 

We note that for T — 0.5 MeV, the minima at lower 
excitation energies are favored. At T = 1.0 MeV, the 
higher energy phase (more broken pairs) is favored. We 
see also, at T = 0.85 MeV, that the free-energy minima 
where we break two and three pairs equal. Where two 



2 



minima coexist, we may have an indication of a phase 
transition. Note however that this is not a phase transi- 
tion in the ordinary thermodynamical sense. There is 
no abrupt transition from a purely paired phase to a 
nonpaired phase. Instead, our system developes several 
such intermediate steps where different numbers of bro- 
ken pairs can coexist. At e.g., T = 0.95 MeV, we find 
again two equal minima. For this case, seniority S = 6 
and S — 8 yield two equal minima. This picture repeats 
itself for higher seniority and higher temperatures. 

If we then focus on the second and third minima, i.e., 
where we break two and three pairs, respectively, the dif- 
ference AF between the minimum and the maximum of 
the free energy, can aid us in distinguishing between a 
first order and a second order phase transition. If AF/N 
remains constant as N increases, we have a second or- 
der transition. An increasing AF/N indicates a first or- 
der phase transition. In Table | we display AF/N for 
N = 10, 12, 14 and 16 at T = 0.85 MeV. It is important 
to note that the features seen in Fig. [I], apply to the cases 
with N = 10, 12 and 14 as well, where T = 0.85 MeV 
is the temperature where the second and third minima 
equal. This means that the temperature where the tran- 
sition is meant to take place remains stable as function 
of number of single-particle levels and particles. This is 
in agreement with the simulations of Lee and Kosterlitz 
p6| . We find a similar result for the minima developed 
at T = 0.95 MeV, where both S = 6 and S = 8 coexist. 
However, due to population inversion, these minima are 
only seen clearly for N — 12, 14 and 16 particles. 

Table Q reveals that AF/N is nearly constant, with 
AF/N ss 0.5 MeV, indicating a transition of second or- 
der. This result is in agreement with what is expected 
for an infinite system. It is also easy to see from Fig. [j], 
that the entropy in the micro canonical ensemble can be 
convex for certain excitation energy ranges, resulting in 
eventual negative heat capacities, as inferred from the 
authors of Refs. [ ^0|JTl| ]. The analysis above however, 
does not lend support to interpreting this as a sign of a 
first order phase transition; see also the recent work of 
Moretto et al. U|. 

We note the important result that for d/G > 1.5, our 
free energy, for N < 16, developes only one minimum for 
all temperatures. That is, for larger single-particle spac- 
ings, there is no sign of a phase transition. This means 
that there is a critical relation between d and G for the 
appearance of a phase transition-like behavior, being a 
reminiscence of the thermodynamical limit. This agrees 
also with e.g., the results for ultrasmall metallic grains 

0. 

We next compute the distributions of zeros (DOZ) of 
the canonical partition function and corresponding poles 
of the specific heat in the complex temperature plane, 
following a scheme outlined by Lee and Yang jl9| , Gross- 
mann et al. |2Cj and Borrmann et al. plj |. Following 
Rcf. pl| , we restrict our discussion to the canonical en- 
semble and denote complex temperatures by B — (3 + it. 
The partition function of Eq. (0) is now a function of B. 



The authors of Ref. M] showed that the distributions 
of zeros are able to reveal the thermodynamic secrets 
of small systems in a distinct manner. Major contribu- 
tions to the specific heat come from zeros close to the 
real axis, and a zero approaching the real axis infinitely 
closely causes a divergence in the specific heat. To char- 
acterize the DOZ close to the real axis, one assumes that 
the zeros lie approximately on a straight line. Defining 
the three parameters r, 7 and a, see Ref. Q for details, 
one can define entirely the nature of the phase transi- 
tion. For a true phase transition in the Ehrcnfcst sense 
we have t — * 0. For this case it has been shown |2C| ] 
that a phase transition is completely classified by a and 
7. In the case a = and 7 = the specific heat CV(/3) 
exhibits a <5-peak corresponding to a phase transition of 
first order. For < a < 1 and 7 = (or 7 7^ 0) the 
transition is of second order. A higher order transition 
occurs for 1 < a and arbitrary 7. 

In Fig. H we show contour plots of the specific heat | 
C V (B) I in the complex temperature plane for N = 11 (a), 
14 (b), and 16 (c) particles at normal pairing d/G = 0.5 
and the N = 14 (d) in the weak pairing limit, d/G = 2. 
The poles are at the center of the dark contour regions. 
We see evidence of two phases in these systems. The 
first phase, labeled I in Fig. 0, is a mixed seniority phase 
while the second phase, //, is a paired phase with zero 
seniority and exists only in even-A systems. No paired 
phase exists in the N = 11 system and no clear bound- 
aries are evident in the weak pairing case. We find that 
for (b) and (c) the DOZ are apparently distributed along 
two lines where the intersection occurs at ri, which is 
the pole closest to the real axis. As the pairing branch 
(for (3 > Pi) only encompasses two points, we are unable 
to precisely determine a along this branch while 7 > 0. 
Based on our free energy results discussed above, we be- 
lieve a along this branch will be positive. In the mixed 
phase branch (for (3 < (3\) we find 7 < 0, and a < in 
all normal-pairing cases. The parameter n, which is a 
measure of discreteness shows a n ~ ^y-i.4±o.i2 Q; e p en _ 
dence. 

We note several significant results of this work. From 
two independent methods we find that the transition 
from the paired seniority zero ground state to a mixed 
phase state is second order. The free-energy analy- 
sis also demonstrates that each transition in seniority 
phases in the microcanonical ensemble is of second or- 
der. The strength of the pairing in these systems de- 
termines the nature of the phase transitions. In par- 
ticular, for a weakly paired system, we found no evi- 
dence for two phases, while normal pairing strengths, 
such as those found in nuclei, may well exhibit the paired- 
phase and mixed seniority phases that we demonstrated 
in this model. We will include more realistic interac- 
tions to investigate this point in future work. We also 
found, using Auxiliary Field Monte Carlo computations 
for this system 22 together with the histogram method 
of Refs. |l6] , |l7[ |, that the energy fluctuations in the canon- 
ical ensemble make it rather difficult to extract useful 
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information on the nature of the phase transitions from 
these techniques. 
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N 


10 


12 


14 


16 


AF/N [MeV] 


0.531 


0.505 


0.501 


0.495 



TABLE I. AF/N for T = 0.85 MeV. See text for further details. 



5 




5 10 15 20 




5 10 15 20 25 



Excitation energy E 



FIG. 1. Free energy from Eq. (§) at T = 0.5, 0.85 and T = 1.0 MeV with d/G = 0.5 with 16 particles in 16 doubly degenerate 
levels. All energies are in units of MeV and an energy bin of 10 -3 MeV has been chosen. 
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FIG. 2. Contour plots of the specific heat in the complex temperature plane for a) N = 11, b) N = 14, and c) N = 16 
particles. Panel d) shows the N — 14 case with weak pairing. The spots indicate the locations of the zeros of the canonical 
partition function. 
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